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Abstract 

Pade approximants are used to find approximate vortex solutions of any wind- 
ing number in the context of Gross-Pitaevskii equation for a uniform condensate and 
condensates with axisymmetric trapping potentials. Rational function and generalised 
rational function approximations of axisymmetric solitary waves of the Gross-Pitaevskii 
equation are obtained in two and three dimensions. These approximations are used to 
establish a new mechanism of vortex nucleation as a result of solitary wave interactions. 
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pulses. 

1 Introduction 

In the last decade the experimental realisation of Bose-Einstein condensation in trapped 
alkali-metal gases at ultralow temperatures has stimulated an intense interest in the produc- 
tion of vortices and vortex arrays and theoretical investigations of their structure, energy, 
dynamics, and stability The condensates of alkali vapours are pure and dilute, so that 
the Gross-Pitaevskii (GP) model which represents the so-called mean-field limit of quan- 
tum field theories gives a precise description of the atomic condensates and their dynamics 
at low temperatures. The same equation has been the subject of extensive studies also in 
the framework of superffuid helium at very low temperature 0, though the high density 
and strong repulsive interactions of superffuid helium restrict the applicability of the GP 
model so that it provides at most a qualitative description. 

In spite of the seeming simplicity of the GP model, which is a defocussing nonlinear 
Schrodinger equation, that has also being extensively studied in other physical systems 
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such as nonlinear optics, not many asymptotic or approximate solutions have been found 
especially for solitary wave solutions in more than one dimensions. Typically one has to 
resort to numerical integration even in the case of a simple straight line vortex. Other 
techniques involve power series expansions that have a small radius of convergence and, 
therefore, are of a very limited use (see more below). The other approach involves quite 
elaborate asymptotic expansions in several different regions and asymptotic matching be- 
tween them ( [W\ |2] ) , which typically yields a useful result for critical parameters of interest 
such as critical velocities, but the resulting approximation of the solution is too compli- 
cated to be used on its own, either in analytic manipulations or as an initial condition 
(possibly with a perturbation) for numerical calculations. Many numerical studies use a 
specific vortex configuration as a starting initial condition (see for instance, |13 |ll4[IT6l l3]) 
and it is desirable to have a simple approximation to the vortex structure that can be used 
in such calculations. 

The reason for the failure of the power series to represent the vortex solution is that 
series diverge in the presence of singularities. There are techniques of a summation theory 
that allow us to overcome this difficulty and represent a given function by a convergent 
expression. In Euler summation this expression is the limit of a convergent series and 
in Borel summation this expression is the limit of a convergent integral. But for these 
methods to work one has to know all the terms of the divergent series exactly before the 
Euler or Borel sum can be found approximately. Pade summation permits us to use only 
a few terms of the divergent series to construct an improved estimate of the function. 

In what follows we shall modify the standard Pade summation method, so that only 
the general forms of the power series at zero and at infinity are used to determine the 
appropriate form of the Pade approximant, but the unknown coefficients will be determined 
recursively from the GP equation. The general idea is that, if the function f(x) has power 
series expansions of the form x n ^2^ =l piX 2l ~ 1 around x = and Qi x ~ 2i a ^ infinity, 

then a Pade approximation of the form 



where = (Zo^iV+n with go = /(°°) will have the same asymptotics at zero and infinity 
as the corresponding power series. 

Our paper is organised as follows. In Section 2 we will derive a Pade approximation 
of the straight line vortex in a uniform condensate for any winding number. Sections 3 
and 4 deal with vortices in a condensate with trapping potentials including an external 
potential with a laser beam. In Sections 5 and 6 we develop the rational and generalised 
rational function approximations of the solitary waves moving with a constant velocity in 
dimensions two and three correspondingly. In Section 7 we study a new mechanism of 
vortex nucleation as a result of solitary wave interactions. 
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2 Straight line vortex in a uniform condensate 



We start with the GP model [SJ 03 in the form 

-2i^ = V^ + (l-|V'| 2 )V', (2) 

where we use dimensionless variables such that the unit of length corresponds to the healing 
length a, the speed of sound is c = l/\/2, and the density at infinity is = 1. The solution 
for the straight line vortex ip = R{r) exp(ira#) with winding number n = 1, 2, ••• in a uniform 
condensate was first obtained by Pitaevskii |17j via numerical integration of the steady GP 
equation 

R"(r) + -R'(r) - ^R{r) + [1 - R 2 {r)}R(r) = 0, (3) 

subject to boundary conditions R(0) = and R(oo) = 1. The asymptotic expansions for 
small r and large r in terms of the power series were obtained by many authors typically 
for n = 1. At small r, the solution can be asymptotically approximated by R(r) ~ 
Y2iLiPi r2i ~ > where the first term of the expansion p\ has to be determined numerically 
(by shooting) as p\ = 0.5827811878 with the rest of the terms then generated recursively as 
p 2 = -0.072847648, p 3 = 0.01128249, p 4 = -0.001781398, • • -. The resulting series is useful 
only within its radius of convergence (i.e., for r < 2.5, 15 ). Asymptotic but divergent 
solution at infinity can be obtained as R(r) ~ S£o^ r_2 *' wnere Qo = 1> Qi = — \, 
q2 = — |, • • •• In view of the expansion in odd powers at zero and even power expansion at 
infinity a Pade approximation of the straight line vortex can be obtained in the form Q 
with n = 1 and = 1 as 

1 + b\r A + 02T 

where we can let 62 = a 2 in the view of the condition p — > 1 as r — > 00. We substitute the 
Pade approximation @ into (jSJ) and expand the resulting expression in series of r setting 
the coefficient at equal powers of r to zero. At C(r 4 ) we get 02 = cii(i>i — 1/4), at C(r 6 ) 
we get 61 = (5 — 32ai)/(48 — 192ai), and at 0(r s ) we get a\ = 11/32 as the positive root 
of 11 + 56ai — 256a^ = 0. The resulting approximation 



, r 2 (0.3437 + 0.0286r 2 ) 
( - r) ~ V 1 + 0.3333r 2 + 0.0286r 4 ^ ^ 

gives the correct asymptotic behaviour at r — ► and at r — > 00 and approximates the 
numerical solution very well everywhere. Figure 1 plots the various approximations to the 
numerically found vortex solution for n = 1 . Note that the procedure described above could 
be formalised by rescaling both equation © and the trial function (@J) by r — > er, equating 
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Figure 1: [Colour online] The plots of the amplitude of the straight vortex line in a uniform condensate 
as a function of distance from the centre of the vortex. The solid black line gives the solution obtained by 
numerically integrating Grey line - the Pade approximation ©. Short dashed and long dashed lines 
give power series expansions at zero (to ©(r 11 )) and at infinity (to 0(r~ 20 )) correspondingly. 
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coefficients at powers of e and setting e = 1 (see also ^5] for power series expansions). 
This procedure will become important when we consider solutions involving more than one 
variable, so that the recursive procedure can be established only by using a different scaling 
of variables. We could also rewrite the equation Q in terms of density p = R 2 (r), so that 
the model equation for the straight line vortex becomes 

d 2 p I dp 1 f dp\ 2 2pn 2 2 

+ rTr ~ YpKlTr) ~^~ 2e {p ~ 1)(> = °' (6) 

A Pade approximation of the straight line vortex with non unit winding number can be 
obtained by observing that R(r) ~ r' n ' at r — > 0, so that we need to consider the expression 

p(r) ~ ^Ka l + a 2{ erf) 

£HoM e r-) 2i + a2(er) 2 M+ 2 

For instance, for \n\ = 2 we obtained 

r 4 (0.0256 + 0.0006264r 2 ) 



(1 + 0.19109r 2 + 0.0196962r 4 + 0.0006264r 6 ) ' 



(8) 
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3 Straight line vortex in a cigar-like trap 

The equation for the density of a condensate in a cigar-like trap with the vortex trapped 
in the centre is 

d 2 p ldp 1 /d/9\2 2p 2 2 2 2 

-7^ + -— -— (— j 2"-2e (p-l + e A r )p = 0, (9) 

ar z r ar 2p\dr / r z 

where A is the dimensionless oscillator frequency. 

For a noninteracting gas, the condensate wave function for a singly quantised vortex 
on the symmetry axis involves the first excited radial harmonic oscillator state V'nonintcr ~ 
r exp(— ^A 2 r 2 ) exp(z#), so we seek a solution of the form 

e 2 r 2 ( Ql +q 2 e 2 r 2 ) 2 2 2 ^ 

PW = : — T 2 2 i A 4 4 eX P(- £ A r )» ( 10 ) 

1 + »ie z r z + foe T - * 

where the Taylor expansion will be taken for the exponential function. The expression (|1U|) 
with 62 = 02 is substituted into © and the terms up to C(e 8 ) are set equal to zero. We 
get 

a 2 = ^ai(-l + 4A 2 + 46i), (11) 

, 5 + 112A 4 - 48A 2 + 32ai (6A 2 - 1) . 
61 " 48 (4 ai + 4A 2 - 1) (12) 

with ai given as a positive root of 

256a 2 (12A 2 - 1) + (4A 2 - 1) 3 (52A 2 - 11) 

- 8ai(384A 6 - 528A 4 + 136A 2 - 7) = 0. (13) 

For example, the solution for A = 0.2 becomes 

0.2833r 2 + 0.0136r 4 A 2n 

p(r) = = jexp(-0.04r 2 . (14) 

Hy 1 l + 0.2581r 2 +0.0136r 4 vy ' y 1 

Similarly to the uniform condensate case, the approximation for the multiply quantised 
vortices can be found in the form 

, , (er) 2 l"l(ai + a 2 e 2 r 2 ) 2 . 2 2l /ir x 

p(r) ~ — rj—^ exp -e X r z . (15) 

El=o^(^) 2i + a2(er) 2 l"l+ 2 
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4 Vortex in an axisymmetric condensate 



A similar procedure can be implemented to find the vortex in a condensate in an axisym- 
metric trap given by the external potential V = X 2 r 2 + X 2 z z 2 . Now, the function p depends 
on two coordinates p = p(r, z), so that the equation for p becomes 

Prr + P -- P ^-% + - 2 ( Pzz -f]-2e\p-l + 6 2 AV + e*\W)p = 0, (16) 
r Ip r z e z \ 2p I 

where we rescaled the variables as r — > er and z — > e 2 z. We seek a solution of the form 
o(r z) - e 2 r 2 ( ai + e 2 a 2 r 2 ) 2 2 2 4 2 2 

We solve the resulting equations to 0(e 8 ), define all the parameters in terms of remain- 
ing two (say ci and a%), that are zeros of coupled polynomials, whose roots can be found 
numerically by fixing A and Xz- For instance, for A = 0.2 and Xz = 0.1 (and setting e = 1) 
we have 

0.232609r 2 + 0.0319011r 4 . n nA 2 n m 2 . , io . 

p = ^ , ^ exp(-0.04r 2 - Q.Qlz 2 ). (18) 

r 1 + 0.342901r 2 + 0.0319011r 4 + 0.00697428z 2 Ly i \ > 

Figure 2 gives the contour plot of this solution in xz— plane. The generalisation to multiply 
charged vortices is carried out in a similar way to previous sections. 

The asymptotics of vortices in a trapped condensate was studied by Konotop and 
Perez-Garcia |12| using a sophisticated multiscale method; their solution was of the form 
Atanh(ar) exp(— br 2 — cz 2 ), where A, a, b, c are some constants that needed to be estimated 
numerically. The use of a Pade approximation not only provides a swift and easy tool for 
generating an approximation of comparable accuracy but also has the ability to include a 
z— dependence that is not confined to the exponential term. 

There are different ways of creating a vortex in a trapped condensate. In particular, 
angular momentum can be transmitted to the condensate by rotationally stirring it with a 
laser beam 0] or by guiding a vortex created at the edge of the condensate by a laser beam 
towards the centre of condensate |2()j . When the vortex is brought to rest at the centre of 
the condensate, the steady external potential can be assumed to have the form 

V{r,z) = X 2 r 2 + X z z 2 + V exp[-r 2 /rf], (19) 

where 77 is the half-width of the laser beam intensity profile. Our procedure for finding 
the approximate solution can be automatically adjusted for finding the vortex density in 
the condensate with potential (J19|). We seek a Pade approximant ()17l) as a solution of the 
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Figure 2: A density contour plot in the xz— plane for a condensate containing a vortex along the z— axis 
as the solution of l|16|l given by l|18^ . The trap parameters are A = 0.2 and Xz ~ 0.1. Luminosity is 
proportional to density, the white area being the most dense. 

20 
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equation 

Ptt l ,-. ■ > ' 1 I ' - - ^ 

- 2e 2 {p - 1 + e 2 AV + e 4 A|z 2 + V exp[-e 2 r 2 /r? ])p = 0. (20) 

In particular, for Vq = ri = 0.8 we get the following solution 

0.18354r 2 + 0.0910069r 4 , , o 

p = ; -. tt exp(-0.04r 2 - O.Olz 2 . (21) 

F 1 + 0.452004r 2 + 0.0910069r 4 + 0.005353z 2 ^ y i \ i 

Figure 3 shows the r— dependent plots of ()18|) and (|21() for various values of z. The form 
of the external potential is given as an inset. Notice how the laser beam causes a slight 
depletion of the condensate close to the centre, followed by an increase in the density 
elsewhere. 



5 Pade approximations of solitary waves in two dimensions 

So far we have shown that Pade approximations are useful for obtaining accurate approxi- 
mations of straight line vortices. In this section we obtain approximate solutions of solitary 
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Figure 3: Plots of the density function p(r, z) at z = 0, 5, 10, 15 for the vortex solutions with f l21H dashed 
line) and without (CHI solid line) the laser beam. The form of the external potential is depicted in the 
inset. 




waves moving with a constant velocity. In a seminal paper, Jones and Roberts |1U| numeri- 
cally determined the entire sequence of solitary wave solutions of the GP equation, such as 
vortex rings, vortex pairs, and rarefaction pulses. Their numerics involved the introduction 
of stretched variables, expansion of the wave function in double Chebyshev-Legendre series 
and a Newton-Raphson iteration of the resulting system of nonlinear algebraic equations. 

We start with finding Pade approximations of solitary waves in two dimensions (2D). 
The pair of two point vortices of opposite circulation centred at (0, yo) and (0, — yo) with 
large yo can be approximated by the superposition of wavefunctions of two straight-line 
vortices of opposite circulation, therefore, by the field with the density 

p = r 2 (Vx 2 + (y- y ) 2 )R 2 (Vx 2 + (y + vo) 2 ), (22) 

and the phase 



S = arctan 


y-yo 


— arctan 


y + yo 




X 




X 



(23) 



where R(r) is given by ©. Figure 4 illustrates how this approximates the solution obtained 
numerically by a Newton-Raphson iteration of the GP equation 

2iU^ = VV + (1 - |Vf)V>, (24) 
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Figure 4: Plots of the density function p(0, y) for vortex pairs with separations 2j/o = 8 (4a) and 2j/o = 1-78 
(4b) taken in a cross-section through the centres of vortices. Solid line - the approximate solution obtained 
by multiplying the wave functions of individual vortices, dashed line - solution obtained numerically. 
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subject to the boundary condition tp — ► 1 as |x| — > oo. Here U is the velocity with which 
the vortices of opposite circulation propel each other in the positive x— direction. 

The approximation (|22|) - (|23|) is sufficiently accurate for large yo, but significantly devi- 
ates from the numerical solution as the separation between the vortices decreases. Another 
difficulty in using this approximation is that one needs to know the separation 2y$ for each 
U in order to construct such an approximation. 

To obtain Pade approximations of 2D solitary solutions we observe that the approxima- 
tion ((2*21) - (|2l?|) written for the real u(x,y) and imaginary v(x,y) parts of the wavefunction 
ip(x,y) = ^/pexp[— iS] = u(x, y) + \v(x, y) are 



u(x,y) = (x 2 + y 2 - yo)R(Vx 2 + (y - y ) 2 )R(^/x 2 + (y + y ) 2 ), 

v(x,y) = -2xy R(^/x 2 + (y - y ) 2 )R(^/x 2 + (y + y ) 2 ), (25) 

where R(r) = R(r)/r. This suggests that a Pade approximation of the solitary wave solution 
moving in x— direction with a constant velocity U can be found in the form 



u(x,y) 



J2 ctijX 2i y 2j 
1 + I] CijX 2i y 2 i 



v ( x , y ) = x y (26) 

1 + 2^ CijX lt y l 3 

We truncate ()26|1 to the lowest order that can give two zeros of the density and unity at 
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Figure 5: The equidistant contour and density plots of the density p — u(x,y) 2 + v(x,y) 2 where u(x,y) 
and v(x, y) are given by (I3U1 . Luminosity is proportional to density, the white area being the most dense. 
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infinity, so we seek an approximation of the form 
u(x,y) = 1 + 
v(x,y) = 



aoo + aiox 2 + aoiy 2 



1 + c 10 x 2 + coiy 2 + c 2 qx a + c u x 2 y 2 + co 2 y 4 
boo + bi x 2 + b iy 2 



1 + ciox 2 + c iy 2 + c 2 qx 4: + c u x 2 y 2 + c 02 y 4 ' 



(27) 



We fix the value of U, substitute (|27|) into the equation (|24|1 . expand in powers of x and 
y and set the coefficients at the first three leading orders to zero. As the result we get 12 
algebraic equations in 11 variables a,ij,bij, and Cij that are compatible and can be solved 
by a computer algebra. The first two leading orders give the analytical expressions of six 
coefficients in terms of remaining five. The sum of squares of the remaining six equations 
is further numerically minimised on the set of remaining five coefficients. 
For U = 0.4 the solution was found as 



1 + 



-1.10048 — 0.095006 ar + 0.01681 y 



1 + 0.3108 x 2 + 0.0192 x A + 0.1030 y 2 + 0.0219 x 2 y 2 + 6.112 x lO" 3 y 4 
x (-0.818647 - 0.06910 x 2 - 0.03756 y 2 ) 



1 + 0.3108 x 2 + 0.0192 x 4 + 0.1030 y 2 + 0.0219 x 2 y 2 + 6.112 x 10- 3 y 4 ' 
The error of the approximation (|28j) is Err\ ~ 0.008, where we defined 

72 ~ 1 M - 2 -v 2 )u) 2 + {-2Uu x + V 2 v + (l-u 2 -v 2 )v) 2 . 



Erri = max(2Uv x + V u + (1 



(28) 



(29) 



From our constructions the expressions (|28|) approximate the solution very well for small x 
and y and give unity for the density at infinity. Two zeros of u(0, y) give the half separation 
between vortices as 0.893997 (compare it with the value found numerically in ^U] as 0.89 
to two significant digits). 
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We can improve the error Err%, by considering not just the global minimum of the 
sum of the squares of 0(e 4 ) error in the power series expansions about the origin, but 
also local minima and choosing the one that minimises Err\. This procedure gives us an 
approximate solution with Err\ < 10~ 4 which is 



-1.09077 - 0.0983212 x 2 - 0.00193044 y 

u = 1 + 



1 + 0.31406 x 2 + 0.02026 x 4 + 0.11677 y 2 + 0.02448 x 2 y 2 + 0.007432 y A 
x (-0.811702 - 0.0717083 x 2 - 0.0465306 y 2 ) 



(30) 



1 + 0.314063 x 2 + 0.02026 x 4 + 0.11677 y 2 + 0.024476 x 2 y 2 + 0.007432 y 4 ' 
The energy, £ , and momentum, p, of solitary wave solutions can be calculated as in 

£ = \ J(l-u 2 -v 2 )(3-2u-u 2 -v 2 )dxdy, 
Up = i /(l -u 2 -v 2 ) 2 dxdy, (31) 



and are £ ~ 8.1 and p ~ 14.2 for u and v given by ()30|) . These values can be compared 
with f Pa 8.16 and p ~ 14.1 found numerically in 10 j. Since we relaxed the precision of the 
approximation for small x and y, the half separation between vortices lost its precision 
as well, becoming 0.87. Figure 5 shows the contour and density plots of the density 
p(x, y) = u(x, y) 2 + v(x, y) 2 of the approximate solution (|3"U|) . 

Finally we consider the third approximation, where we further relax the precision of 
the approximation at the origin, but impose additional constraints at large x and y. The 
asymptotic solution for large x and y was obtained in JOJ as 



. . m(U - \m)x 2 - mU{l - 2U 2 )y 2 
u(x, y) f« 1 H 



(x 2 + (1 - 2U 2 )y 2 ) 2 
v(x,y) « - x 2 + ™ 2U * )y 2 > W-°°> (32) 

where m is a constant termed 'the stretched dipole moment' of the wave because of the 
factor 1 — 2U 2 . If we would like our solution to have this asymptotic behaviour at infinity, 
we have to let 

aio = C2om(U — |m), aoi = — C2om(l — 2U 2 )U, bio = —mc2o 

boi = -mc 2 o(l-2U 2 ), c n = 2c 20 (l - 2U 2 ), c 02 = c 20 (l - 2U 2 ) 2 . (33) 

The expansion of the equation (|24jl about zero to C(e 2 ) gives four equations on seven un- 
knowns U, m, aoo, boo, C20> cio, coi- We can solve these equations analytically for c 2 o, cio, cqi, 
and boo- In particular, the following approximate relation between the minimum of the real 
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Figure 6: Equidistant contour plot of the density p = u(x,y) 2 + v(x,y) 2 where u(x,y) and v(x, j/)are given 
by J37J. 

Figure was not included in this submission to keep the size below the maximum allowed. 

part of the wavefunction u(0, 0) = aoo + 1 and the slope of the imaginary part at the origin 
in the direction of a solitary motion v x (0,0) = boo was found 

„ l( 0,0) = W0,0)-l)±5±^^. (34) 

Next we fix one parameter, say U, to specify one particular solution of the family, and find 
aoo an d m that minimise the integral error 

Err 2 = j (2Uv x + V 2 u + (1 - u 2 - v 2 )u) 2 + (-2Uu x + V 2 v + (1 - u 2 - v 2 )v) 2 dV, (35) 

where dV = dxdy. 

By implementing this procedure, we obtained the following approximation of the vortex 
pair moving with U = 0.4 

. , , -1.14026 - 0. 150112 x 2 -0.0294564 y 2 

v ' yj 1 + 0.35022 x 2 + 0.03032 x 4 + 0.15905 y 2 + 0.04123 x 2 y 2 + 0.01402 y 4 ' 

x (-0.830953 - 0.108296 x 2 -0.073641 y 2 ) 
v ( x,y > 1 + 0.35022 x 2 + 0.03032 x 4 + 0.15905 y 2 + 0.04123 x 2 y 2 + 0.01402 y 4 ^ ' 

with the stretched dipole moment found as m = 3.57 (compared with the value 3.55 found 
numerically in jlOj). The absolute error of this approximation is Err\ ~ 0.0018. 

In order to find the solitary wave with a single zero of the wave function, we implement 
the latter procedure by fixing aoo = — 1- This will guarantee that the intersection of zeros 
of real and imaginary parts is only at the origin. We found that U = 0.45, m = 3.32 and 
the approximate solution with Err\ = 0.00089 is 

. , , -1 - 0.137233 x 2 -0.0303092 y 2 

v ; l + 0.36838x 2 + 0.03397x 4 + 0.15803y 2 + 0.04082x 2 y 2 + 0.01227 y 4 

, . x (-0.803361 -0.112913 x 2 - 0.0678506 y 2 ) , 

v (x v) = - — — (37) 

v ' 1 + 0.36838x 2 + 0.03397 x 4 + 0.15803 y 2 + 0.04082 x 2 y 2 + 0.01227 y 4 ' v ; 

On Figure 6 the contour plot of p(x, y) is given. This solution is interesting because it 
represents the borderline case between vortex pair and vortex-free solutions. 
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A Pade approximation of a rarefaction pulse can be obtained similarly. For instance, 
if U = 0.5 we get aoo = —0.826 and m = 3.1, so that the solution is 

. . , -0.825937 -0.114393a; 2 -0.0271467 y 2 

nix v) = 1 ~\~ 

v ; 1 + 0.37355 x 2 + 0.03495 x 4 + 0.14235 y 2 + 0.03495 x 2 y 2 + 0.008737 y 4 ' 

, s x (-0.737901 - 0.108587 x 2 -0.0542934 y 2 ) , 

v (x v) = - — - d8) 

v ,y; 1 + 0.37355 x 2 +0.03495 x 4 +0.14235 y 2 + 0.03495 x 2 y 2 + 0.008737 y 4 ' v ; 

We will use this approximation of the rarefaction pulse in Section 7, where we study vortex 
nucleation. 



6 Generalised rational function approximation of axisym- 
metric solitary waves in three dimensions 

To determine approximations of the axisymmetric vortex rings and rarefaction pulses mov- 
ing along the x— axis with the constant velocity U in three dimensions (3D), we need to 
solve (|24[) with the Laplacian written in cylindrical coordinates 

- 2Uv x = u xx + u ss + \u s + (1 - u 2 - v 2 )u, 
2Uu x = v xx + v ss + ±v s + (1 - u 2 - v 2 )v , (39) 



where s = \Jy 2 + z 2 . The asymptotics of the solitary waves at large distances from the 
origin was found in ^0] to have the form 

2mUx 2 - 2mU(l - 2U 2 )s 2 

U[X,S ' ~ 1 + (X 2 + (1 - 2C/2) s 2)5/2 ' 

<*>*) * - {x 2 + { i-2u 2 )s 2 r^ |x| ^°°- (40) 

It is clear that approximation by rational functions can not behave at large distances as 
(|40j) . therefore, we will use 'generalised' rational functions. In particular, we determined 
that among many possibilities the following expressions give a reasonable approximation 
everywhere 

/ x 1 , aoo + aio^ 2 + a is 2 + mc^XJ(2x 2 - (1 - 2XJ 2 )s 2 ) 
u(x,y) = H 



(1 + c 10 x 2 + cois 2 + c 20 (x 2 + (1 - 2U 2 )s 2 ) 2 y/ i ' 

6oo + b 10 x 2 + b 01 y 2 - m^\x 2 + (1 - 2U 2 )s 2 ) 2 
1 ,!J> (1 + c l0 x 2 + c 0l s 2 + c 20 (x 2 + (1 - 2U 2 )s 2 ) 2 )V^ 1 ) 

From the two leading order expansions of (|39j) about the origin we determine a±o, aoi, box, b±o 
and the rest of the unknowns will be found by minimising (|35j) . where dV = sdsdx, for 
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a fixed U. Among many local minima that give small Err 2^ we choose the one with a 
stretched dipole moment close to the one obtained in JU]. In particular, the small vortex 
ring with U = 0.6 and Err\ = 0.003 was found as ()41j) with 

a 00 = -1.1, a i = 0.0170524, a 10 = 0.0289452, m = 8.97 

600 = -0.953, 6 i = -0.0049767, b w = -0.0594346, (42) 

c i = 0.04, cio = 0.21, c 20 = 0.00612. 

We calculate the energy and momentum of the solitary wave solutions in 3D as (see |10j ) 

p = 2tt J ((u — l)v x — vu x )s ds dx, (43) 

to get £ ~ 58.8 and p ~ 78 for (|41|1 with (|43|) . that can be compared with the values found 
numerically in jlOj : £ ~ 56.4 and p ~ 78.9. The radius of the ring given by (|41|) - ()43f) is 
1.059 (with 1.06 found numerically). 

Finally we give our approximation of the rarefaction pulse moving with U = 0.63 found 
on the lower branch of the dispersion curve calculated numerically in jlOj : 

a 00 = -0.79792, a 01 = 0.00388059, a w = 0.00882276, m = 8.37 
6 00 = -0.7981, 6 i = -0.012783, b w = -0.0574092, 

c i = 0.0399, cio = 0.199, c 20 = 0.0058. (44) 

For this approximation of the rarefaction pulse (|41 |l -(|44j l we have £ ~ 54.4, p ~ 72.2, 
and m = 8.37 (compared with numerical values found in [TJ]]: £ ~ 52.3, p ~ 72.2, and 
m = 8.37). 



7 Vortex nucleation 

The approximations developed in the previous sections allow one to study interactions 
among the solitary wave solutions of the GP equation. Without an accurate starting point 
in numerical calculations it would be impossible to separate clearly the effect of interactions 
from the evolution of each solitary wave by itself as it settled down from a poor initial guess. 

In this section we will use the approximations developed above to show that vortex pairs 
and vortex rings can appear as a result of an interaction among the solitary wave solutions 
of the GP equation. Previously, the nucleation of vortices in a uniform condensate has 
been connected to critical velocities |2] or instabilities of the initial states TJ. We will 
show that interactions between various, even vortex-free, solitary waves result in energy 
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Figure 7: [Colour online] The snapshots of the contour plots of the density cross-section of a condensate 
obtained by numerically integrating the GP model (3- Initial condition is ip(t = 0) = ty(x — 5, 0)' i l l (x + 5, 0), 
where "if — u + iv, with u and v given by 1411 - 144(1 . Black solid lines show zeros of real and imaginary 
parts of tp, therefore, their intersection shows the position of topological zeros. Both low and high density 
regions are shown in darker shades to emphasise intermediate density regions. Only a portion of an actual 
computational box is shown. 
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and momentum transfer that can lead to vortex nucleation. Rarefaction pulses on the 
lower branch of the dispersion curve have lower energy and momentum than vortex rings, 
therefore, such rarefaction pulse may evolve into a vortex ring if interactions with other 
solutions add enough energy and momentum to the rarefaction pulse. 

This scenario is supported by direct numerical simulations, performed with the same 
numerical method as in our previous work [21 El- In these computations we follow the 
evolution of two rarefaction pulses moving in the computational box of dimensions Z) 3 = 
80 3 . The faces of the box are open to allow sound waves to escape; this is achieved 
numerically by applying the Raymond-Kuo technique |18j . 

We prepare our initial conditions by superimposing the wavefunctions ipi of solitary 



15 



Figure 8: [Colour online] The isosurface pi poo = 0.15 of a condensate. Initial condition is ip(t = 0) = 
^(x — 5, 0)^(a; + 5, 0), where 9 — u + iv, with u and w given by l)41[l - H44|p . At t = the solution is vortex 
free. AT t = 22.5 the front pulse contains a closed topological zero. At t = 40 a well-formed vortex ring 
has developed with the back solution becoming too shallow to be shown for this isosurface. Only a portion 
of an actual computational box is shown. 




wave solutions of the GP equation found in the previous sections, tp(t = 0) = Y\ ipi. If the 
distance between such solutions is large enough, then such a superposition will not lead to 
a significant initial sound emission and solitary waves will preserve their form initially. Our 
first initial state consists of two rarefaction pulses (|41j) - (|44[) positioned the distance 10 apart 
and moving towards each other. Two pulses collide and pass through each other without 
loss of energy. Next we will create a field nonuniformity by placing two rarefaction pulses a 
distance 10 apart that move in the same direction. This time the effect two solitary waves 
have on each other is non-symmetric. As a result, the rarefaction pulse moving behind 
transfers part of its energy and momentum to the pulse moving at front, so that the latter 
transforms into a vortex ring and slows down, whereas the former spreads out and speeds 
up. This process leads to an even closer interaction of the two solitary waves and an even 
more rapid transfer of energy from the solitary wave that moves behind to the one moving 
at front. Eventually almost all of the energy and momentum of the former is transferred 
to the latter, which becomes a vortex ring of energy and momentum that are only slightly 
less than twice the energy and momentum of each of the initial rarefaction pulses. The 
remaining small energy is emitted as sound waves. Figure 7 gives the graphical illustration 
of this process through the snapshots of the density cross-sections. Figure 8 shows the 
density isoplots at |^| 2 = 0.15 of the various stages of the ring formation. 

Similarly, energy and momentum transfer takes place between different types of solitary 
waves. In our next calculation we start with a rarefaction pulse followed by a large vortex 
ring that moves in the same direction. Initially both the distance between these two solitary 
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waves and the radius of the ring were taken to be 10. The rarefaction pulse is moving faster 
than the vortex ring, so the distance between them rapidly increases. Nevertheless, there is 
an energy and momentum transfer that allows the rarefaction pulse to evolve into a vortex 
ring of a small radius and the large vortex ring to shrink slightly. Apart from a small loss 
of energy to sound waves, the total energy of these two solitary waves is almost conserved 
throughout this transformation. These processes are shown in graphical form through the 
contour plots of the density cross-sections (Figure 9) and the density isoplots at \ip\ 2 = 0.1. 
(Figure 10). 

Figure 9: [Colour online] The snapshots of the contour plots of the density cross-section of a condensate 
obtained by numerically integrating the GP model (5J. The initial condition is ij)(t = 0) = "ifi(x, 0)^>2(x — 
5, s+ 10)^2 (a! - 5, s — 10), where i&i = u+iv, with u and v given by JIT)- (3D and * 2 = R(\/x 2 +y 2 )e ie with 
R given by @. Black solid lines show zeros of real and imaginary parts of ip, therefore, their intersection 
shows the position of topological zeros. Both low and high density regions are shown in darker shades to 
emphasise intermediate density regions. Only a portion of an actual computational box is shown. 




Similar transfer of energy and momentum from one solitary wave to another takes place 
in 2D. In particular, there is a transfer of energy between solitary waves moving in the 
same direction, whereas two colliding solitary waves interact elastically. Nevertheless, this 
elasticity of interactions can be broken by introducing other solitary waves. As an example, 
we show the evolution of an initial condition consisting of two colliding rarefaction pulses 
in close vicinity of a widely separated vortex pair. As a result of the interaction, another 
vortex pair is created and the resulting two pairs of vortices move apart in direction making 
small angles with the positive x— axis. Figure 11 shows the snapshots of the density of a 
condensate at various moments of time as the solution evolves. 
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Figure 10: [Colour online] The isosurface pf poo = 0.1 of a condensate. The initial condition is ip(t = 
0) = *i(a;,0)*2(a; — 5, s + 10) — 5, s — 10), where *i = u + iv, with u and v given by 101}- GU an d 
* 2 = R(Vx 2 + s 2 )e w with i? given by ©. At f = the front pulse is vortex free. AT t — 22.5 the front 
pulse contains a closed topological zero. At t = 40 the front pulse evolved into a well-formed vortex ring. 
Only a portion of an actual computational box is shown. 
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Figure 11: [Colour online] The snapshots of the contour plots of the density of a condensate obtained 
by numerically integrating the CP model (|5J in 2D. Initial condition is if>(t = 0) = ipi(x — 10, 0)^1 (£ + 
10, 0)ip2(x,y + 4)4>2( x ,y — 4)) where tpi — u + iv, with u and v given by 1381 and tp2 = R(\/ x 2 + y 2 )e lB with 
R given by JSj- Black solid lines show zeros of real and imaginary parts of ip, therefore their intersection 
shows the position of topological zeros. Both low and high density regions are shown in darker shades to 
emphasise intermediate density regions. Only a portion of an actual computational box is shown. 

Figure was not included in electronic submission to keep the size below maximum 

allowed. 

8 Conclusions 

In summary we have presented a new technique for finding approximate vortex solutions 
of the GP equation in a uniform condensate and in condensates with axisymmetric traps. 
These solutions have simple analytic expressions, correct asymptotic behaviours at zero and 
infinity and approximate the entire solutions quite well elsewhere. We envision that the 
use of such approximations will allow one to set up accurate initial vortex configurations 
for numerical calculations and will make explicit analytic manipulations possible. 

We have also developed a technique for obtaining approximations of the solitary solu- 
tions such as vortex pairs and rarefaction pulses in two and three dimensions. The found 
approximations are shown to give a very low error and have simple analytical form. These 
approximations are used to elucidate the energy and momentum transfer between different 
solitary wave solutions by direct numerical integration of the GP equation. The process of 
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vortex nucleation is one of the consequences of such a transfer. 
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